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ABSTRACT 

How  to  reduce  the  computational  error  is  a  key  issue  in  numerical  modeling  and  simulation.  The  higher  the 
order  of  the  difference  scheme,  the  less  the  truncation  error  and  the  more  complicated  the  computation.  For 
compromise,  a  simple,  three-point  accuracy  progressive  (AP)  finite-difference  scheme  for  numerical  calculation 
is  proposed.  The  major  features  of  the  AP  scheme  are  three-point,  high-order  accuracy,  and  accuracy  progressive. 
The  lower-order  scheme  acts  as  a  “source”  term  in  the  higher-order  scheme.  This  treatment  keeps  three-point 
schemes  with  high  accuracy.  The  analytical  error  estimation  shows  the  sixth-order  accuracy  that  the  AP  scheme 
can  reach.  The  Fourier  analysis  of  errors  indicates  the  accuracy  improvement  from  lower-order  to  higher-order 
AP  schemes.  The  Princeton  Ocean  Model  (POM)  implemented  for  the  Japan/East  Sea  (JES)  is  used  to  evaluate 
the  AP  scheme.  Consider  a  horizontally  homogeneous  and  stably  stratified  JES  with  realistic  topography.  Without 
any  forcing,  initially  motionless  ocean  will  keep  motionless  forever;  that  is  to  say,  there  is  a  known  solution 
(V  =  0).  Any  nonzero  model  velocity  can  be  treated  as  an  error.  The  stability  and  accuracy  are  compared  with 
those  of  the  second-order  scheme  in  a  series  of  calculations  of  unforced  flow  in  the  JES.  The  three-point  sixth- 
order  AP  scheme  is  shown  to  have  error  reductions  by  factors  of  10-20  compared  to  the  second-order  difference 
scheme.  Due  to  their  three-point  grid  structure,  the  AP  schemes  can  be  easily  applied  to  current  ocean  and 
atmospheric  models. 


1.  Introduction 


Most  ocean  models  are  hydrostatically  balanced.  Such 
a  balance  sets  up  a  lowest  limit,  (Ax,  A y)L,  for  the  grid 
spacings.  If  the  model  resolution  is  finer  than  that  limit, 
nonhydrostatic  ocean  models  should  be  used.  If  the  hy¬ 
drostatic  balance  is  kept  unchanged,  the  grid  spacing 
cannot  be  artificially  small  and  the  reduction  of  dis¬ 
cretization  errors  within  the  grid  spacing  limitation  be¬ 
comes  important.  A  logical  approach  is  to  use  high-order 
difference  schemes. 

High-order  schemes  can  be  applied  to  the  cr-coordi- 
nate  ocean  models  for  error  reduction.  Let  (x*,  y*,  z) 
denote  Cartesian  coordinates  and  (x,  y,  cr)  cr  coordinates. 
The  conventional  relationships  between  the  two  coor¬ 
dinates  are 

a*  =  x,  y*  =  y,  z  =  crH(x,  y), 


where  H  is  the  water  depth,  and  z  and  cr  increase  ver¬ 
tically  upward  such  that  z  =  cr  =  0  at  the  surface  and 
cr  =  —  1  and  z  =  ~H  at  the  bottom.  The  horizontal 
pressure  gradient  becomes  a  difference  between  two 
terms  using  the  cr  coordinates, 


dp  dp  cr  dH  dp 
Sx*  dx  H  dx  da 


(1) 


Corresponding  author  address:  Dr.  Peter  Chu,  Dept,  of  Ocean¬ 
ography,  Naval  Postgraduate  School,  Monterey,  CA  93943. 

E-mail:  chu@nps.navy.mil 


which  leads  to  a  so-called  hydrostatic  inconsistency  with 
a  large  truncation  error  at  a  steep  topography  (Gary 
1973;  Haney  1991).  Mellor  et  al.  (1994)  argue  that  if 
the  pressure  gradient  force  is  computed  as  a  vertical 
integral  of  discretized  Jacobian  of  density  and  z  coor¬ 
dinate,  the  problem  of  hydrostatic  inconsistency  can  be 
avoided  and  the  spurious  residual  pressure  gradient  is 
drastically  reduced.  Furthermore,  Song  (1998)  shows 
that  a  modified  Jacobian  formulation  on  the  discredited 
pressure  gradient  term  with  a  specially  designed  non¬ 
linear  weighting  procedure  drastically  reduces  the  res¬ 
idue  errors. 

Thus,  we  have  two  approaches  to  reduce  errors  of  the 
cr-coordinate  ocean  models:  1)  bringing  certain  sym¬ 
metries  of  the  continuous  forms  into  the  discrete  level 
to  ensure  cancellations  of  these  terms  (e.g.,  Mellor  et 
al.  1994;  Song  1998;  Song  and  Wright  1998),  and  2) 
increasing  numerical  accuracy  (e.g.,  McCalpin  1994; 
Chu  and  Fan  1997,  1998,  1999,  2000).  Note  that  the 
use  of  high-order  schemes  does  not  conflict  with  the 
use  of  symmetric  form.  We  may  combine  the  two  ap¬ 
proaches  (symmetric  forms  with  high-order  schemes)  to 
further  diminish  the  computational  errors. 

Usually,  the  higher  the  accuracy  of  the  scheme,  the 
more  the  grid  points  are  needed.  For  example,  a  C-grid 
fourth-order  ordinary  scheme  needs  four  grid  points 
(McCalpin  1994)  and  a  C-grid  sixth-order  ordinary 
scheme  needs  six  grid  points  (Chu  and  Fan  1997)  for 
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computing  the  first  derivative.  However,  the  current 
ocean  models  generally  use  three-point  second-order 
schemes.  Conversion  of  the  existing  ocean  models  from 
second-  to  high-order  (fourth  or  sixth  order)  schemes 
becomes  easier  if  the  high-order  scheme  preserves  the 
three-point  grid. 

Concurrently,  Hirsh  (1975)  and  Kreiss  (1975)  have 
proposed  Hermitian  compact  techniques  using  fewer 
grid  points  (three  instead  of  five)  at  each  grid  point  to 
solve  partial  differential  equations  (PDEs).  Later  on,  as 
pointed  out  by  Adam  (1977),  the  truncation  errors  are 
usually  four  to  six  times  smaller  than  the  same  order 
noncompact  schemes.  Since  then,  much  work  has  been 
done  in  developing  compact  schemes  for  various  ap¬ 
plications,  such  as  an  implicit  compact  fourth-order  al¬ 
gorithm  (Navon  and  Riphagen  1979),  a  fourth-order 
compact  difference  scheme  for  nonuniform  grids  (Goe- 
dheer  and  Potters  1985),  and  fourth-order  and  sixth- 
order  compact  difference  schemes  for  the  stagged  grid 
(Chang  and  Shirer  1985).  Recently,  Chu  and  Fan  (1998, 
1999,  2000)  proposed  a  three-point  sixth-order  com¬ 
bined  compact  difference  (CCD)  scheme  for  uniform, 
nonuniform,  and  staggered  grids.  This  scheme  follows 
earlier  work  on  use  of  second  derivatives  in  compact 
differences  (such  as  Rubin  and  Khosla  1977), 

i  (aji+k  +  bj;+k  +  ckflk)  =  0, 

k=  —  l 

which  is  referred  as  the  Hermitian  formula.  Here,  /  is 
a  dependent  variable.  The  CCD  schemes  use  implicit 
solvers  to  calculate  various  differences  (twin-tridiagonal 
solver)  and  to  solve  PDEs  (triple-tridiagonal  solver). 
Such  a  complicated  procedure  may  limit  the  application 
of  the  CCD  schemes  to  the  ocean  models. 

In  this  study,  we  propose  a  simple,  three-point  ac¬ 
curacy  progressive  (AP)  finite-difference  scheme  for  nu¬ 
merical  calculation.  The  major  features  of  the  AP 
scheme  are  three  point,  high-order  accuracy,  and  ac¬ 
curacy  progressive.  The  lower-order  scheme  acts  as  a 
“source”  term  in  the  higher-order  scheme.  This  treat¬ 
ment  keeps  in  three-point  grid  with  high  accuracy  and 
is  explicit  in  time  integration. 


2.  Progress  of  accuracy 
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the  second-order  differences  of  /'  and  f"  at  the  grid 
point  x„ 
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are  the  source  terms  of  the  standard  Pade  (fourth  order) 
scheme  in  the  right-hand  sides  of  (1)  and  (2). 

Following  this  path,  we  propose  a  series  of  three- 
point  AP  finite-difference  schemes  with  /',  /",  .  .  .  ,  f(t> 
computed  by  a  lower-order  scheme  as  source  terms  in 
higher-order  schemes.  They  have  the  following  features: 
three-point,  accuracy  progressive,  and  only  one  tridi¬ 
agonal  solver  needed.  The  second-order  AP  scheme  is 
the  same  as  the  central  difference  scheme,  and  the 
fourth-order  AP  scheme  is  the  same  as  the  standard  Pade 
scheme. 


3.  AP  schemes 

a.  Standard  grid 

Let  (Af/Ax,  A2f/Ax2),  ( 8f/8x ,  82f/8x 2),  and  ( Df/Dx , 
D2flDx2)  be  the  second-order,  fourth-order,  and  sixth- 
order  differences  for  (/',  /"),  respectively.  The  accuracy 
progress  from  the  second-order  scheme  into  the  fourth- 
order  AP  scheme  is  given  by 


Let  the  dependent  variable  /(x)  be  defined  on  the 
interval,  0  <  x  s  L.  Use  a  uniform  grid,  0  =  x,  <  x2 
<  x3  <  ■  •  •  <  xN  <  xN+l  =  L  with  a  spacing  h  = 
xi+1  —  x,  =  UN.  Let  the  dependent  variable  /(x)  at  any 
grid  point  x,  and  two  neighboring  points  x,  ,  and  x,+1 
be  given  by  /;,  , ,  and  fi+1  and  let  its  derivatives  at 

the  two  neighboring  points  x:  ,  and  xi+1  be  given  by 

fi-u ru  •  •  • ,  m,  and /;+1, r;+1, . . . ,  /«. 

A  new  approach  is  to  relate  higher-order  schemes  to 
lower-order  schemes.  For  example,  in  the  standard  Pade 
(fourth  order)  scheme. 


(5) 


+  a2 


+ 


(6) 


Taylor  series  expansion  shows  that  if  a1  =  1/4,  /S,  = 
3/2,  a2  =  1/10,  and  /32  =  6/5,  the  schemes  (5)  and  (6) 
have  the  fourth-order  accuracy. 

The  accuracy  progress  from  the  fourth-order  AP 
scheme  into  the  sixth-order  AP  scheme  can  be  defined  by 
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When  a,  =  -7/254,  /3,  =  120/127,  y,  =  -17/254,  a2 
=  —5/94, 1 62  =  —102/47,  and  y2  =  144/47,  the  schemes 
(13)  and  (14)  have  the  sixth-order  accuracy. 
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j82  =  3,  and  y2  =  —9/4,  the  schemes  (7)  and  (8)  have 
the  sixth-order  accuracy  (see  appendix). 


b.  Staggered  C  grid 

The  second-order  AP  scheme  on  a  staggered  C  grid 
is  given  by 


A/\  1 

/  =r(L-/,-), 

Ax/i+i/2  h 


A  2f\  1 

A i).  -  - 2f' +  /-)- 
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The  fourth-order  scheme  is  given  by 
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Taylor  series  expansion  shows  that  if  a,  =  1/22,  /3,  = 
12/11,  a2  =  1/10,  and  /32  =  6/5,  the  schemes  (11)  and 
(12)  have  the  fourth-order  accuracy. 

The  sixth-order  AP  scheme  on  a  staggered  C  grid  can 
be  defined  by 
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4.  Fourier  analysis  of  errors 

Fourier  analysis  of  errors  is  commonly  used  to  eval¬ 
uate  various  difference  schemes,  described  extensively 
in  Orzag  (1971),  Kreiss  and  Oliger  (1972),  Lele  (1992), 
and  Chu  and  Fan  (1998).  As  pointed  out  by  Lele  (1992), 
Fourier  analysis  provides  an  effective  way  to  quantify 
the  resolution  characteristics  of  differencing  approxi¬ 
mations. 

For  the  purpose  of  Fourier  analysis,  the  dependent 
variable  f{x)  is  assumed  to  be  periodic  over  the  domain 
[0,  L\  of  the  independent  variable,  that  is,  fl  =  fN+l 
and  h  =  L/N.  The  dependent  variable  may  be  decom¬ 
posed  into  Fourier  series, 

k=NI2 

fix)  =  2  fke^'L\  (15) 

k=—(N/2) 

where  i  =  V  — 1.  It  is  convenient  to  introduce  a  scaled 
wavenumber  w  =  2TrkhlL  =  lirk/N,  and  a  scaled  co¬ 
ordinate  s  =  xlh.  The  Fourier  modes  in  terms  of  these 
are  simply  exp(/ws).  The  exact  first-order  and  second- 
order  derivatives  of  (15)  generate  a  function  with  exact 
Fourier  coefficients: 


fl 


However,  the  Fourier  coefficients  of  the  derivatives  ob¬ 
tained  from  the  differencing  scheme  might  not  be  the 
same  as  the  exact  Fourier  coefficients,  that  is. 


(fl)fd  =  ifDfd  = 

where  w'  =  w'{w)  and  w"  =  w"{w)  are  the  modified 
wavenumber  (both  real  numbers)  for  the  first-order  and 
second-order  differencing.  The  smaller  the  difference 
between  the  exact  and  modified  wavenumbers,  the  better 
the  difference  scheme. 

According  to  Lele  (1992),  if  the  modified  wave- 
numbers  of  the  current  generalized  difference  schemes  as 

/;  + «(/;+:  +  fid  =  afi+i  ~  fi-' 


ri  + «(/ 
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the  modified  wavenumbers  will  be 
a  sinw 

w  (w)  =  -  and 

1  +  2  a  cosw 
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(a)  First  derivative  approximation  (b)  Second  derivative  approximation 


Polar  plot  of  phase  speed  anisotropy  for  the  first  derivative  approximation; 
wavenumber(magnitude):  1/50,  5/50, .  45/50,  50/50 


Fig.  1.  Fourier  analysis  of  error  for  derivative  approximation:  (a)  second-order  AP  scheme 
(second-order  central  scheme),  (b)  fourth-order  AP  scheme  (standard  Pade  scheme),  (c)  sixth- 
order  combined  compact  scheme,  (d)  sixth-order  AP  scheme,  and  (e)  exact  differentiation. 


For  the  second-order  AP  scheme  (central  difference 
scheme),  we  have 

Wj  (w)  =  sinw,  w"  (w)  =  [2(1  —  cosw)]1/2,  (18) 


for  the  fourth-order  AP  scheme  (standard  Pade’ 
scheme),  we  have 


w[{w) 


3  sinw 
2  +  cosw’ 


w"(w) 


12(1  —  cosw) 

5  +  cosw 

09) 


and  for  the  sixth-order  AP  schemes  (7)  and  (8),  we  have 
63  +  27  cosw 

w6(w)  = 


(8  +  7  cosw)(5  +  cosw) 


sinw 


(20) 


wl  (w)  = 


4[  6  —  6  cosw  —  ~We(w)  sinw 


COSW 


(21) 


Among  various  difference  schemes,  the  modified 
wavenumbers  of  the  first-order  differencing  w'  (Fig.  la) 
and  of  the  second-order  differencing  w"  (Fig.  lb)  of  the 
AP  scheme  are  closest  to  the  exact  wavenumber  w.  Be¬ 
sides,  we  computed  the  errors  of  the  modified  waven¬ 
umbers  for  the  first-order  derivative  (Fig.  lc) 

lw'(w)  —  w  I 

error  =  - 

w 

and  for  the  second-order  derivative  (Fig.  Id) 

I  w"(w)  —  wl 

error  =  - , 

w 

respectively.  We  clearly  see  the  accuracy  increase  of  the 
AP  schemes  from  the  second-order  to  the  fourth-order, 
and  from  the  fourth-order  to  the  sixth-order.  Further¬ 
more,  the  performance  of  the  sixth-order  AP  scheme  is 
very  close  to  the  CCD  scheme. 
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(a) 


1  2  3  4  N-2N-1  N  N+1 

I _ I _ I _ I _  _____  _ ! _ I _ ! _ I 

J-* -  Variable  - *| 

|* -  2nd  AP  - *| 

|< -  4nd  AP  - *| 

U -  6nd  AP  - J 


(b) 

N-2N-1  N  N+1 

    »  l  «  l  »  l  »  l 

Variable  - 

2nd  AP  - *j 

4nd  AP  - >| 

6nd  AP  - >| 

Fig.  2.  Boundary  conditions  for  AP  schemes:  (a)  standard  grid  and 
(b)  staggered  C  grid. 


In  multidimensional  problems  the  phase  error  of  first- 
order  differencing  scheme  appear  in  the  form  of  an¬ 
isotropy: 


(C'p)fd(w,  9)  = 


w'(w,  9) 
w 


(cos0)w'(w  cos0)  +  (sin0)v/(tv  sin0) 
w 

(22) 


Figure  le  shows  polar  plots  of  phase  speed  anisotropy 
of  various  schemes  for  first  derivative  approximations. 
The  phase  speeds  for  wavenumbers  (magnitude)  whr  = 
1/50,  5/50,  .  .  . ,  45/50,  50/50  are  plotted.  Here,  we  also 
see  that  the  sixth-order  AP  scheme  shows  improvement. 


5.  Boundary  treatment 

Due  to  the  progressive  accuracy,  the  AP  schemes  need 
special  treatment  near  the  boundaries  x,  and  xN+1.  For 
the  standard  AP  schemes,  the  derivatives  can  be  com¬ 
puted  at  (x2,  x3, .  .  . ,  Xjy)  using  the  second-order  scheme, 
at  (x3,  x4,  .  .  . ,  xN  , )  using  the  fourth-order  scheme,  and 
at  (x4,  x5,  .  .  .  ,  xN_2)  using  the  sixth-order  scheme  (Fig. 
2a).  Thus,  we  use  the  second-order  scheme  for  the  grid 
points  x2  and  xN  (next  to  the  boundaries),  the  fourth- 
order  scheme  for  the  grid  points  x,  and  xN  , ,  and  the 
sixth-order  scheme  for  the  grid  points  (x4,  x5, .  .  . ,  xN_2). 
For  the  sixth-order  AP  scheme  on  the  staggered  C  grid, 
we  have  the  similar  treatment  (Fig.  2b). 


132  134  136  138  140  142 

Longitude  (E) 

Fig.  3.  Geography  and  isobaths  showing  the  bottom  topography  of 
the  JES  model. 
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realistic  ocean  model,  the  analytical  solution  is  hard  to 
find.  Consider  a  horizontally  homogeneous  and  stably 
stratified  ocean  with  realistic  topography.  Without  forc¬ 
ing,  initially  motionless  ocean  will  keep  motionless  for¬ 
ever;  that  is  to  say,  we  have  a  known  solution  (V  =  0). 
Any  nonzero  model  velocity  can  be  treated  as  an  error. 

We  use  the  Princeton  Ocean  Model  (POM;  Blumburg 
and  Mellor  1987),  implemented  for  the  Japan/East  Sea 
(JES),  to  evaluate  the  AP  scheme.  The  model  contains 
the  treatment  (symmetric  forms  for  term  cancellation) 
proposed  by  Mellor  et  al.  (1994). 


b.  JES  geography 

The  Japan  Sea,  known  as  the  East  Sea  in  Korea,  has  a 
steep  bottom  topography  that  makes  it  a  unique  semi- 
enclosed  ocean  basin  overlaid  by  a  pronounced  monsoon 
surface  wind.  The  JES  covers  an  area  of  106  km2.  It  has 
a  maximum  depth  in  excess  of  3700  m,  and  is  isolated 
from  open  oceans  except  for  small  (narrow  and  shallow) 
straits,  which  connect  the  JES  with  the  North  Pacific 
through  the  Tsushima/Korean  and  Tsugaru  Straits  and  with 
the  Okhotsk  Sea  through  the  Soya  and  Tatar  Straits.  The 
smoothed  JES  bathymetry  is  shown  in  Fig.  3. 


6.  Evaluation  using  the  Princeton  Ocean  Model 

a.  Known  solution 

Accuracy  of  any  numerical  scheme  should  be  eval¬ 
uated  by  either  an  analytical  or  known  solution.  For  a 


c.  Model  description 

The  POM  is  a  time-dependent,  cr-coordinate,  primi¬ 
tive  equation  circulation  model  on  a  three-dimensional 
grid  that  includes  realistic  topography  and  a  free  surface 
(Blumburg  and  Mellor  1987).  The  model  contains  94 
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Longitude  (E) 

Fig.  4.  Horizontal  pressure  gradient  force  error  I A I  (unit:  N)  at  day  20,  using  (a)  the  second- 
order  scheme  without  the  symmetric  cancellation,  (b)  the  second-order  scheme  with  the  symmetric 
cancellation,  and  (c)  the  sixth-order  AP  scheme  with  the  symmetric  cancellation. 


X  100  X  15  horizontally  fixed  grid  points.  The  hori¬ 
zontal  spacing  of  10'  latitude  and  longitude  (approxi¬ 
mately  1 1.54-15.18  km  in  the  zonal  direction  and  18.53 
km  in  the  latitudinal  direction)  and  15  vertical  cr-co- 
ordinate  levels.  The  model  domain  is  from  35.0°  to 


51.5°N,  and  from  127.0°  to  142.5°E.  The  bottom  to¬ 
pography  is  obtained  from  the  smoothed  Naval  Ocean¬ 
ographic  Office  Digital  Bathymetry  Data  Base  5  min 
by  5  min  resolution.  The  horizontal  diffusivities  are 
modeled  using  the  Smagorinsky  (1963)  form  with  the 
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(a)  (b) 


Average  IAI  (N)  Maximum  IAI  (N) 


Fig.  5.  Vertical  dependence  of  (a)  horizontally  averaged  IAI  and 
(b)  maximum  I A I  (unit:  N)  at  day  20  using  the  second-order  scheme 
without  the  symmetric  cancellation  (a:  dash-dotted),  the  second-order 
scheme  with  the  symmetric  cancellation  (b:  solid),  and  the  sixth-order 
AP  scheme  with  the  symmetric  cancellation  (c:  dashed). 


coefficient  chosen  to  be  0.2  for  this  application.  No 
atmospheric  forcing  is  applied  to  the  model. 

Closed  lateral  boundaries,  that  is,  the  modeled  ocean 
bordered  by  land,  were  defined  using  a  free-slip  con¬ 
dition  for  velocity  and  a  zero-gradient  condition  for  tem¬ 
perature  and  salinity.  No  advective  or  diffusive  heat, 
salt,  or  velocity  fluxes  occur  through  these  boundaries. 
At  open  boundaries,  we  use  the  radiative  boundary  con¬ 
dition  with  zero  volume  transport. 

d.  Initial  conditions 

The  model  was  integrated  with  all  three  components 
of  velocity  (u,  v,  w)  initially  set  to  zero,  and  with  the 


(a) 


same  potential  density  used  by  Mellor  et  al.  (1994)  for 
investigating  the  horizontal  pressure  gradient  error, 

f  -  1  028  -  0.003  exp(^),  (23, 

which  can  be  treated  as  an  approximation  to  the  area- 
averaged,  vertical  density  distribution.  Here  z  is  the  ver¬ 
tical  coordinate,  Hp  =  1000  m,  and  p0  =  1 0 ;  kg  m~3. 

e.  Mode  splitting 

For  computational  efficiency,  the  mode  splitting  tech¬ 
nique  (Blumburg  and  Mellor  1987)  is  applied  with  a 
barotropic  time  step  of  25  s,  based  on  the  Courant- 
Friederichs-Levy  (CFL;  1928)  computational  stability 
condition  and  the  external  wave  speed;  and  a  baroclinic 
time  step  of  900  s,  based  on  the  CFL  condition  and  the 
internal  wave  speed. 

f.  Symmetric  cancellation 

The  POM  contains  the  symmetric  cancellation 
scheme  (Mellor  et  al.  1994).  Without  the  cancellation 
scheme,  we  compute  the  scalar  pressure  field  first  by 
vertical  integration  of  density  and  then  its  horizontal 
gradient  via  straightforward  discretization.  With  the 
cancellation  scheme,  we  compute  the  pressure  gradient 
force  as  the  vertical  integral  of  discretized  Jacobian  of 
density  and  z  coordinate.  With  this  treatment,  the  prob¬ 
lem  of  hydrostatic  inconsistency  can  be  avoided  and  the 
spurious  residual  gradient  is  drastically  reduced  (Mellor 
et  al.  1994). 


(b) 


Fig.  6.  Volume  transport  streamfunction  (m3  s  ')  using  the  second-order  scheme  on  (a)  day  5  and  (b)  day  20. 
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(a) 


(b) 


Fig.  7.  Volume  transport  streamfunction  (m3  s  ')  using  the  sixth-order  AP  scheme  on  (a)  day  5  and  (b)  day  20. 


7.  Horizontal  pressure  gradient  force  error 

When  the  density  field  (23)  is  used,  the  horizontal 
pressure  gradient  should  be  zero,  theoretically;  however, 
it  is  usually  nonzero  in  the  cr-coordinate  system  due  to 
the  computational  error.  We  may  use  the  horizontal  pres¬ 
sure  gradient  force  on  the  three-dimensional  grid  cell 

A  =  AxAy(HAcr)(dp/dx dp/dyf) 

at  day  20  to  evaluate  the  performance  of  the  three 
schemes  in  computing  the  horizontal  pressure  gradient 
vector  for  the  JES  POM:  the  second-order  scheme  with¬ 
out  the  symmetric  cancellation  (Fig.  4a),  the  second-order 
scheme  with  the  symmetric  cancellation  (Fig.  4b),  and 
the  sixth-order  AP  scheme  with  the  symmetric  cancel¬ 
lation  (Fig.  4c)  at  four  different  cr  levels  (1,  5,  9,  and 
14).  Here,  Ajc,  Ay,  and  Act  are  the  grid  spacings  in  the 
x,  y,  and  cr  axes.  The  a  levels  1  and  14  are  the  surface 
and  bottom,  respectively.  We  use  the  cubic-spline  inter¬ 
polation  to  compute  dp/da.  Using  the  same  second-order 
difference  scheme,  we  find  an  order  of  magnitude  re¬ 
duction  in  the  horizontal  pressure  gradient  error  with  the 
symmetric  cancellation  (Fig.  4b)  compared  to  without 
the  symmetric  cancellation  (Fig.  4a).  Using  the  same 
symmetric  cancellation  (Mellor  et  al.  1994),  we  find  an¬ 
other  order  of  magnitude  reduction  in  the  horizontal  pres¬ 
sure  gradient  error  with  the  sixth-order  AP  scheme  (Fig. 
4c)  compared  to  the  second-order  scheme  (Fig.  4b). 

We  compute  horizontally  averaged  and  maximum  val¬ 
ues  of  I A I ,  over  a  cr  level,  to  compare  the  performance 
of  these  schemes  (Fig.  5).  The  symmetric  cancellation 
scheme  (Mellor  et  al.  1994)  has  the  capability  to  dras¬ 
tically  increase  the  accuracy  using  the  same  second- 
order  scheme.  For  example,  the  horizontally  averaged 
(maximum)  values  of  I A I  generally  reduce  by  factors 


of  5-32  (8-23)  with  the  symmetric  cancellation  scheme 
compared  to  without  the  symmetric  cancellation 
scheme.  Moreover,  using  the  symmetric  cancellation  we 
see  a  further  reduction  in  the  horizontally  averaged 
(maximum)  values  of  I A I  by  factors  of  6-9  (4-9)  with 
the  sixth-order  AP  scheme  compared  to  the  second-or¬ 
der  scheme. 

8.  Error  estimation  using  the  JES  POM 

We  integrated  the  diagnostic  JES  POM  with  Mellor 
et  al.’s  (1994)  symmetric  cancellation  scheme  from  zero 
velocity  and  horizontally  integrated  density  field  (23) 
for  20  days  with  the  second-order  scheme  and  the  sixth- 
order  AP  scheme,  respectively.  Without  computational 
error,  the  model  ocean  should  keep  motionless.  With  the 
computational  error,  the  model  generates  false  motion. 
The  absolute  value  of  the  velocity  is  taken  as  the  error 
velocity.  The  smaller  the  error  velocity,  the  more  ac¬ 
curate  the  numerical  scheme  would  be. 

a.  Volume  transport  streamfunction 

Volume  transport  streamfunction  (1P)  on  days  5  and 
20  are  presented  in  Fig.  6  for  the  second-order  scheme 
and  in  Fig.  7  for  the  sixth-order  AP  scheme.  All  the 
figures  show  a  multi-eddy  structure  in  coincidence  with 
the  steep  bottom  topography  (Fig.  3).  The  difference 
between  two  \P  isolines  indicates  the  volume  transport 
at  that  location.  The  smaller  the  difference,  the  less  the 
volume  transport.  The  contour  interval  of  'P  isolines  is 
five  times  smaller  in  Fig.  7  than  in  Fig.  6.  This  indicates 
the  error  reduction  by  a  factor  of  5  using  the  AP  sixth- 
order  scheme  compared  to  using  the  second-order 
scheme. 
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Time  (day) 


Fig.  8.  Temporally  varying  (a)  peak  and  (b)  global  error  velocities  for 
the  second-order  (solid)  and  the  sixth-order  AP  (dashed)  schemes. 


b.  Temporal  variations  of  global  and  peak  error 

velocities 

Two  quantities  are  commonly  used  for  identifying  the 
error  evolution:  peak  and  global  error  velocities.  The 
peak  error  velocity  is  the  maximum  absolute  value  the 
spurious  velocity  generated  by  the  pressure  gradient  er¬ 
rors.  Figure  8a  shows  the  time  evolution  of  the  peak 
error  velocity  for  the  first  20  days  of  integration  with 
the  second-order  and  sixth-order  AP  schemes.  The  peak 
error  velocities  increase  to  a  maximum  value  and  then 
decay  as  the  inertial  oscillation  superimposed  into  some 
mean  values:  4.80  X  10~3  m  s_1  (the  second-order 
scheme),  and  0.62  X  10~3  m  s_1  (the  sixth-order  AP 
scheme).  The  maximum  peak  error  velocity  is  5.38  X 
10~3  m  s_I  (0.76  X  10~3  m  s_1)  at  day  12  (day  7.5) 
for  the  second-order  (sixth-order  AP)  scheme.  The 
steady  approach  of  the  peak  error  velocities  to  these 
values  for  the  three  schemes  indicates  the  stability  of 
the  computation.  Furthermore,  the  peak  error  velocity 


for  the  second-order  case  is  roughly  eight  times  that  of 
the  sixth-order  AP  case. 

At  each  time  step,  we  compute  the  spatially  averaged 
error  velocity  over  the  whole  JES  (called  the  global  error 
velocity)  to  represent  the  overall  accuracy  of  the 
schemes.  Figure  8b  shows  the  time  evolution  of  the 
global  error  velocity  for  the  first  20  days  of  integration 
with  the  second-order  and  the  sixth-order  AP  schemes. 
It  increases  with  time  at  a  much  slower  pace  using  the 
sixth-order  AP  scheme  than  the  second-order  scheme. 
On  day  20  the  global  error  velocity  is  2.02  X  10~4  m 
s-1  for  the  second-order  scheme  and  0.34  X  10~4  m  s_1 
for  the  sixth-order  AP  scheme.  The  global  error  velocity 
for  the  second-order  scheme  is  roughly  six  times  that 
of  the  sixth-order  AP  scheme. 

c.  CPU  time 

We  integrate  the  POM  for  20  days  on  the  SGI  Origin- 
200  machine  to  identify  the  cost  of  pressure  gradient 
schemes.  The  central  processing  unit  (CPU)  times  are 
21.067  h  for  the  second-order  scheme  and  21.167  h  for 
the  sixth-order  AP  scheme,  respectively.  There  is  no 
noticeable  CPU  time  increase  from  the  second-order 
scheme  to  the  sixth-order  AP  scheme. 

9.  Open  boundary  influence 

a.  Coastal  geometry 

One  of  the  difficult  problems  in  shallow  water  mod¬ 
eling  is  the  uncertainty  of  the  open  boundary  condition 
(OBC).  At  open  boundaries  where  the  numerical  grid 
ends,  the  fluid  motion  should  be  unrestricted.  Ideal  open 
boundaries  are  transparent  to  motions.  Two  approaches, 
local  type  and  inverse  type,  are  available  for  determining 
OBC  (Chu  et  al.  1997)  with  low-order  conditions.  Be¬ 
fore  converting  any  ocean  model  from  a  second-order 
scheme  to  a  sixth-order  scheme,  it  is  important  to  verify 
if  a  high-order  interior  plus  low-order  boundary  con¬ 
ditions  (such  as  open  boundary  conditions)  would  de¬ 
grade  the  interior  solution.  Consider  a  horizontally  ho¬ 
mogeneous  and  stably  stratified  coastal  ocean  with  a 
longitudinal  and  straight  coastline  and  three  open 
boundaries  in  the  south  (20°N),  the  north  (33.5°N),  and 
the  east  (x  =  xE  =  320  km).  Choose  coordinates  so  that 
the  y  axis  coincides  with  the  coast,  positive  x  pointing 
offshore.  The  coastal  water  is  divided  into  three  parts 
(Fig.  9a):  shelf  (0  <  x  <  x0),  slope  (x0  <  x  <  x,  = 
200  km),  and  deep  water  (xt  <  x  ^  xE ). 

Here,  the  water  depth  ( h )  is  only  a  function  of  offshore 
distance  x.  At  the  coast  (x  =  0),  the  water  has  a  minimum 
depth  h„tm  (Fig.  9).  An  analytical  bottom  topography  is 
proposed  such  that  shelf  and  slope  are  arcs  of  two  cir¬ 
cles.  The  shelf  has  a  smaller  radius  (r),  and  the  slope 
has  a  larger  radius  ( R ).  The  two  arcs  are  connected  in 
such  a  way  that  the  tangent  of  the  bottom  topography, 
dh/dx,  is  continuous  at  the  shelf  break  (x  =  x0).  This 
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Fig.  9.  Coastal  geometry  with  open  boundaries:  (a)  horizontal  view,  (b)  cross-coastal  view,  and  (c) 

three-dimensional  view. 


requirement  is  met  using  the  same  expanding  angle  (6) 
for  both  arcs  (Fig.  9b).  Thus,  this  bottom  topography 
has  5  degrees  of  freedom:  hmin,  /?max,  r,  R,  and  ft 
Off-shore  distance  x  and  bottom  topography  h  are 
nondimensionalized  by 

sinf?  x  li(x)  —  h  ■ 

= _ k*  =  — — _ — 

1  -  cos  dXy  hmax  -  hmin 

and  in  turn  r  and  R  are  nondimensionalized  by 

_  _ 1 _ 

~~  (1  +  k)(l  -  cos ey 


(1  +  k)(l  -  COS0)'  r' 


The  bottom  topography  (Fig.  9)  can  be  represented 
analytically  by 

h*(x*) 

r *  —  (r*2  —  x*2)1/2,  if  x*  ^  x* , 

1  -  R* 

+  [/?*2  —  (x*  —  xf)2]1/2,  if  x^  <  x*  <  xf, 

1 ,  if  x*  >  xf . 

As  in  the  previous  section,  initially  motionless  coastal 
ocean  with  an  exponentially  stratified  density  (23)  and 
without  forcing  will  keep  motionless  forever,  that  is,  V 
=  0.  Any  nonzero  model  velocity  can  be  treated  as  an 
error.  We  use  the  POM  implemented  for  such  a  coastal 
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(a)  (b) 


Fig.  10.  Volume  transport  streamfunction  (m3  s  ')  on  day  10  using  (a)  the  second-order  scheme  and  (b) 

the  sixth-order  AP  scheme. 


region  to  test  the  open  boundary  influence  for  the  high- 
order  AP  scheme.  At  open  boundaries,  we  use  the  ra¬ 
diative  boundary  condition  with  zero  volume  transport. 
The  parameters  are  chosen  in  this  study  by 

0  =  60°,  k  =  5,  h,mn  =  1  m,  /zmax  =  2000  m. 

This  implies  that  the  width  of  the  shelf  break  (x0)  is 
33.33  km. 

In  the  numerical  integration,  the  space  is  discretized  by 

Ax  =  4  km.  Ay  =  20  km.  Act  =  1/14, 

and  barotropic  and  baroclinic  time  steps  are  chosen  as 
6  and  180  s  in  order  to  satisfy  the  CFL  condition. 

b.  Volume  transport  streamfunction 

Volume  transport  streamfunction  ( "'T/ )  on  day  10  is  pre¬ 
sented  in  Fig.  10a  for  the  second-order  scheme  and  in  Fig. 
10b  for  the  sixth-order  AP  scheme.  Both  figures  show  a 
multi-eddy  structure  in  coincidence  with  the  continental 
slope  (Fig.  9).  The  difference  between  two  4'  isolines 
indicates  the  volume  transport  at  that  location.  The  smaller 
the  difference,  the  less  the  volume  transport.  The  contour 
interval  of  4'  isolines  is  around  200  times  smaller  in  Fig. 
10b  than  in  Fig.  10a.  This  indicates  the  error  reduction  by 
a  factor  of  200  using  the  AP  sixth-order  scheme  compared 
to  using  the  second-order  scheme. 

c.  Horizontal  pressure  gradient  force  error 

We  compute  horizontally  averaged  and  maximum  val¬ 
ues  of  IAI  at  day  10  over  a  a  level,  to  compare  the 
performance  of  the  sixth-order  AP  and  second-order 


schemes  (Fig.  11).  The  horizontally  averaged  (solid)  and 
maximum  (dashed)  I A I  values  are  quite  small  at  the 
upper  ocean  for  both  the  second-order  scheme  (Fig.  11a) 
and  the  sixth-order  AP  scheme  (Fig.  lib).  Using  the 
second-order  scheme,  the  errors  (IAI  values)  increase 
monotonically  versus  cr,  with  largest  errors  next  to  the 
ocean  bottom  (cr  =  14):  2943  N  for  the  averaged  IAI 
and  18  582  N  for  the  maximum  IAI.  Using  the  sixth- 
order  AP  scheme  (Fig.  1  lb),  the  errors  (IAI  values)  are 
much  smaller  than  the  second-order  scheme.  The  largest 
errors  occur  at  one  level  higher  than  the  second-order 
scheme  (cr  =  13):  44  N  for  the  averaged  IAI  and  977 
N  for  the  maximum  IAI. 

Error  ratios  between  second-order  and  sixth-order 
schemes  show  drastic  reductions  by  factors  of  30-514 
in  the  horizontally  averaged  I A I  values  (solid)  and  by 
factors  of  43-1160  in  the  horizontally  maximum  IAI 
values  (dashed)  (Fig.  11c).  The  error  reduction  becomes 
very  evident  in  the  middle  level  (cr  =  6)  of  the  ocean. 

d.  CPU  time 

We  integrate  the  POM  for  this  open  boundary  test 
case  for  10  days  on  the  SGI  Origin-200  machine  to 
identify  the  cost  of  pressure  gradient  schemes.  The  CPU 
times  are  501.07  min  for  the  second-order  scheme  and 
504.19  min  for  the  sixth-order  AP  scheme,  respectively. 
There  is  no  noticeable  CPU  time  change  from  the  sec¬ 
ond-order  scheme  to  the  sixth-order  AP  scheme. 

10.  Conclusions 

1)  This  study  shows  a  process  of  building  a  series  of 
accuracy  progress  difference  schemes  with  the  lower- 
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(a)  2nd-ord  scheme 


(b)  6th-ord  ap  scheme 


(c)  2nd-ord  over  6th-ord  ap  scheme 


5  10  15 

IAI  (103  N) 


Fig.  11.  Vertical  dependence  of  mean  (solid)  and  maximum  (dashed)  horizontal  pressure  gradient  force 
errors  IAI  (unit:  N)  on  day  10  (a)  using  second-order  scheme  and  (b)  using  sixth-order  AP  scheme,  (c) 
Error  ratio  between  second-order  scheme  vs  sixth-order  AP  scheme.  The  symmetric  cancellation  is  used  in 
both  cases. 


order  scheme  as  the  “source”  term  in  the  higher-order 
scheme.  Such  schemes  have  the  following  features: 
three-point,  accuracy  progressive,  and  only  tridiagonal 
solver  needed.  The  second-order  AP  scheme  is  the  same 
as  the  central  difference  scheme,  the  fourth-order  AP 
scheme  is  the  same  as  the  standard  Pade  scheme,  and 
the  sixth-order  AP  scheme  has  the  sixth-order  accuracy. 

2)  Fourier  analysis  shows  the  progressive  error  de¬ 
crease  in  modified  wavenumber  from  the  second-order 
AP  scheme  to  the  fourth-order  scheme,  and  then  from 
the  fourth-order  scheme  to  the  sixth-order  scheme. 

3)  The  Princeton  Ocean  Model  implemented  for  the 
Japan/East  Sea  is  used  to  demonstrate  the  benefit  of 
using  the  sixth-order  AP  scheme.  Two  calculations  of 
unforced  flow  are  performed  with  the  second-order 
scheme  and  the  sixth-order  AP  scheme,  respectively. 
The  results  show  that  the  sixth-order  AP  scheme  has 
error  reductions  by  factors  of  5-8  compared  to  the  sec¬ 
ond-order  difference  scheme. 

4)  Using  the  sixth-order  AP  scheme  for  computing 
the  pressure  gradient  does  not  require  much  more  CPU 
time.  Taking  POM  for  the  Japan/East  Sea  as  an  example, 
the  CPU  time  for  the  sixth-order  AP  scheme  is  almost 
the  same  as  for  the  second-order  scheme. 

5)  Both  standard  and  staggered  AP  schemes  were  de¬ 
veloped.  The  sixth-order  AP  scheme  on  a  staggered  C  grid 
can  be  easily  implemented  into  current  ocean  models. 

6)  The  high-order  AP  scheme  has  the  capability  to 
handle  the  open  boundary  condition.  There  is  no  deg¬ 
radation  of  the  interior  solution  using  high-order  AP 
schemes  and  low-order  boundary  conditions. 

7)  Future  studies  include  applying  AP  schemes  to 
nonuniform  grid  systems  as  well  as  designing  even  high¬ 
er-order  schemes  such  as  an  eighth-order  AP  scheme. 
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APPENDIX 


Error  Estimation  of  AP  Schemes  on  a 
Standard  Grid 

The  truncation  error  is  usually  obtained  from  Taylor 
series  expansion.  The  first  and  second  derivatives  at  the 
grid  point  x ,  are  represented  by 


/'  =  ^(/,+l  -  /, ,)  +  om  =  +  om, 

/,"  =  ^(/;+l  -  2/,  +  /,.,)  +  om 


+  om. 


which  shows  the  second-order  accuracy  of  A/7Ax,  A2// 
Ax2.  Also,  /'  and  /"  at  the  grid  point  x,  can  be  rep¬ 
resented  by  (standard  Pade  scheme) 


/,'  +  1  +  /'- 1) 
/,"  +  i  +  fl  i) 


-  — — +  om, 

2  2  h 

6  fj+i  -  2/,  + 

5  h2 


f(6) 

+  —  h4  +  om, 

200 


which  indicates  the  fourth-order  accuracy;  that  is,  the 
difference  between  the  first-order  derivative  and  its 
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fourth-order  difference  counterpart  at  two  grid  points 
xi+ 1  and  x,-  ,  should  be 


-  |1J.)  =0(/t4)  and  (Al) 


=  0(/t4).  (A2) 

Subtraction  of  (A2)  from  (Al)  leads  to 


* f 

8x 


=  h(fi+l  +  +  o(h y 


The  difference  between  the  second-order  derivative  and 
its  fourth-order  difference  at  two  grid  points  x,+1  and 
x j _ [  should  be 


-  (  —  )  =  — /z4  +  <9(/z6)  and  (A3) 

'Sx2/.,,  200 


/'Li  -  (— )  =  —  h4  +  <9(/z6).  (A4) 

J‘  1  \8x1)i_l  200  '  v  ' 

Subtraction  of  (A4)  from  (A3)  leads  to 


sy 

Sx2 


sy 

Sx2 


=  M/'Ll  -  /'-,)  +  0(/*6). 

(A5) 


Chu  and  Fan  (1998)  provide  the  following  error  es¬ 
timation: 


/;  +  ^(/;+ 1  +  -  ygtCi  -  /,/.) 

=  y^'+1  ^  ^  j  -  0(Jib)  and  (A6) 

^(/:+,  -  /; ,)  -  g  +  /," 

=  1  ~  +  ^i+1)  -  °(/l 6)-  (A7) 

Elimination  of  (f"+1  —  /"  ,)  from  (A5)  and  (A6)  leads  to 


/;  +  ^(/;+1  -  /;  i) 

=  i5//,+1  +  /,  A  h_\(sy\  _  !sy\ 

8  \  2 h  )  16[^x2j;+1  [Sx2J:  I 

-  0(hb).  (A8) 

Equations  (A8)  and  (A7)  indicate  sixth-order  accuracy 
for  the  AP  schemes  (7)  and  (8). 
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